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The 2-d spin i Heisenberg antiferromagnet with exchange coupling J is investigated on a periodic 
square lattice of spacing a at very small temperatures using the loop-cluster algorithm. Monte Carlo 
data for the staggered and uniform susceptibilities are compared with analytic results obtained in 
the systematic low-energy effective field theory for the staggered magnetization order parameter. 
The low-energy parameters of the effective theory, i.e. the staggered magnetization density M s — 
0.30743(l)/a 2 , the spin stiffness p s = 0.18081(11) J, and the spin wave velocity c = 1.6586(3)Jo 
are determined with very high precision. Our study may serve as a test case for the comparison of 
lattice QCD Monte Carlo data with analytic predictions of the chiral effective theory for pions and 
nucleons, which is vital for the quantitative understanding of the strong interaction at low energies. 
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Introduction — Partly motivated by the relation of 
antifcrromagnetism to high-temperature superconductiv- 
ity, during the past twenty years quantum spin models, 
such as the spin i Heisenberg antiferromagnet on the 
square lattice, have been studied in great detail. Since 
this system is strongly coupled, numerical simulations 
play an important role in its quantitative analysis. In 
this way, it has been shown that the SU(2) S spin symme- 
try breaks down spontaneously to its U(l) s subgroup at 
zero temperature. As a result, massless Goldstone bosons 
— the antiferromagnetic magnons — dominate the low- 
energy physics. The magnon dynamics can be described 
quantitatively using a low-energy effective field theory for 
the staggered magnetization order parameter [lHll • Low- 
energy phenomena can then be investigated analytically, 
order by order in a systematic derivative expansion. 

Systematic effective field theories also play an impor- 
tant role in the low-energy physics of the strong interac- 
tion. On the one hand, lattice QCD describes the un- 
derlying dynamics of quarks and gluons beyond pertur- 
bation theory, but can only be investigated by very large 
scale Monte Carlo calculations. On the other hand, chiral 
perturbation theory, the systematic low-energy effective 
field theory for pions — the pseudo- Goldstone bosons of 
the spontaneously broken SU(2)l X SU(2)r chiral sym- 
metry of QCD — has been investigated analytically in 
great detail. The predictions of the effective theory de- 
pend on a number of low-energy parameters, including 
the pion decay constant, the chiral condensate, as well 
as the higher-order Gasser-Lcutwylcr coefficients. For 
the quantitative understanding of the strong interaction 
at low energies, it is of central importance to accurately 
determine the values of the low-energy parameters by 
comparison of lattice QCD Monte Carlo data with an- 
alytic chiral perturbation theory predictions. In recent 
years, there has been substantial progress in this direc- 
tion, and the leading-order low-energy parameters have 



been determined with a few percent accuracy. Extending 
this to the higher-order low-energy parameters, as well as 
reaching higher precision while keeping complete control 
of systematic errors will be a major challenge for lattice 
QCD in the near future. 

The 2-d spin i Heisenberg antiferromagnet can serve 
as an ideal test case, in which the interplay between high- 
precision numerical simulations of the underlying micro- 
scopic system and high-order calculations in the cor- 
responding systematic low-energy effective field theory 
can be investigated quantitatively. In contrast to lattice 
QCD which is much more complicated, the Heisenberg 
model can be simulated with very efficient methods, and 
has been inve stig ated in several high-accuracy numeri- 
cal studies [oHlo] The first very precise determination 
of the low-energy constants of the 2-d spin | Heisen- 
berg antiferromagnet was performed in [(| using the loop- 
cluster algorithm This study was based on a cubical 
space-time geometry, for which the inverse temperature 
p = l/T, which determines the extent of Euclidean time, 
is compatible with the spatial size L, i.e. (3c « L. By com- 
parison of the Monte Carlo data with analytic 2-loop re- 
sults of Hasenfratz and Niedermayer, obtained in the sys- 
tematic low-energy effective field theory for the staggered 
magnetization order parameter [f| , the staggered magne- 
tization density M s , the spin stiffness p Sl and the spin 
wave velocity c were determined as M. s — 0.3074(4)/a 2 , 
p s = 0.186(4)J, and c = 1.68(l)Ja. A few years later, 
the development of the continuous-time simulation tech- 
nique Q enabled numerical investigations of quantum 
spin models at very low temperatures. This allowed a 
comparison of Monte Carlo data with analytic 1-loop 
results in the cylindrical space-time regime at very low 
temperatures (3c > L, which led to M s = 0.3083(2)/a 2 , 
p s = 0.185(2)J, and c = 1.68(l)Ja, in statistical agree- 
ment with the results obtained in the cubical space-time 
regime. The fit in the cylindrical regime required an 
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incorporation of 2-loop corrections with adjustable pre- 
factors, because these effects had not been determined 
analytically at that time. Recently, Niedermayer and 
Weiermann have closed this gap by performing the corre- 
sponding analytic 2-loop calculation in the effective the- 
ory [l8j|. Slab-like space-time geometries with L > /3c 
have been investigated in Q. In that study, using finite- 
size scaling, very long spatial correlation lengths up to 
350000 lattice spacings have been investigated. A com- 
bined fit of Monte Carlo data in the cubical, cylindri- 
cal, and slab geometries then gave M s = 0.30797(3)/a 2 , 
p s = 0.1800(5)J, and c = 1.657(2) Ja. In a recent 
study using a zero-temperature valence-bond projector 
method, Sandvik and Evertz obtained the very accurate 
result M s = 0.30743(l)/a 2 . Although the discrepancy 
between these two results for M s is at the per mille 
level, it is statistically significant. Indeed, in a high- 
precision analysis of the constraint effective potential of 
the staggered magnetization, which relied on 2-loop pre- 
dictions in the effective theory by Gockeler and Leutwyler 
1|| 2(|, we suspected that the previously obtained esti- 
mates M s = 0.3083(2)/a 2 Q and M s = 0.30797(3)/a 2 
@, which were dominated by Monte Carlo data in the 
cylindrical regime, are afflicted by an underestimated 
systematic error resulting from a truncation of the See- 
ley expansion described in (BJ. In this paper, we re- 
turn to the cylindrical regime and clarify the discrep- 
ancy. This will result in a confirmation of the value 
M s = 0.30743(l)/a 2 obtained in as well as in a de- 
termination of p s = 0.18081(11)J and c = 1.6586(3)Ja 
with unprecedented precision. In the cubical regime we 
determine c by tuning /3 until the squares of the spatial 
and temporal winding numbers become identical. In this 
way, the uncertainties of both M s and p s resulting from 
the fits are drastically reduced. 

While reaching fractions of a per mille precision for 
the low-energy parameters may seem unnecessary from a 
condensed matter physics perspective, it is reassuring for 
the ongoing efforts to combine lattice QCD with chiral 
perturbation theory in order to accurately determine 
the fundamental low-energy parameters of the strong 
interaction. Using the 2-d Heisenberg model as a test 
case, our analysis demonstrates that very precise Monte 
Carlo data combined with 2-loop effective field theory 
predictions for a variety of physical quantities indeed 
leads to a completely consistent very high precision de- 
termination of the fundamental low-energy parameters. 

Microscopic Model and Corresponding Observ- 
ables — The spin | Heisenberg model considered in this 
study is defined by the Hamilton operator 
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where 1 and 2 refer to the two spatial unit- vectors. Fur- 
ther, J in cq.([T]) is the antifcrromagnctic exchange cou- 



pling. A physical quantity of central interest is the stag- 
gered susceptibility 
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dt iTr[M s 3 (0)M s 3 (t)exp(-/3 J ff)] 
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Here Z = Trexp(— j3H) is the canonical partition func- 
tion. The staggered magnetization order parameter is de- 
fined as M s = J2x(~ l) Xl+X2 S x . Another relevant quan- 
tity is the uniform susceptibility 
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dt -Tr[M 3 (0)M 3 (£)exp(-/3£f)]. (3) 



Here M = S x is the uniform magnetization. Both 
Xs and Xu can be measured very accurately with the 
loop-cluster algorithm using improved estimators [(§■ In 
particular, in the multi-cluster version of the algorithm 
the staggered susceptibility is given in terms of the 
cluster sizes \C\ as Xs — jj^ (£c l^l 2 )- Similarly, the 

uniform susceptibility X u = & (W?) = (£ c W t {C) 2 ) 
is given in terms of the temporal winding number 
Wt = £ c Wt(C) which is the sum of winding numbers 
Wt(C) of the loop-clusters C around the Euclidean time 
direction. Similarly, the spatial winding numbers are 
defined by W t = £ c W ^ C ) with i e i 1 * 2 }- 

Low-Energy Effective Theory for Magnons - 

Due to the spontaneous breaking of the SU(2) S spin 
symmetry down to its U(l) s subgroup, the low-energy 
physics of an antiferromagnet is governed by two mass- 
less Goldstone bosons, the magnons. A systematic low- 
energy effective field theory for magnons was developed 
in jlj-jj]. The staggered magnetization of an antifer- 
romagnet is described by a unit-vector field e(x) that 
takes values in the coset space SU(2) s /U(l) s = S 2 , 
i.e. e(x) = (ei(x) , e2{x) , es(x)) with e(x) 2 = 1. Here 
x = (x±,X2,t) denotes a point in (2+l)-dimensional 
space-time. To leading order, the Euclidean magnon low- 
energy effective action takes the form 



S[e] = / dx\ I dx2 

Jo Jo Jo 



(It 



x y \ d\e- die + d 2 e- d 2 e+ ^d t e- <9 t ej , (4) 

where t refers to the Euclidean time-direction. It should 
be noted that the effective field theory described by cq. ((4]) 
is valid as long as the conditions Log ^> 1 and (icp s 3> 1 
are satisfied. As demonstrated in ja], once these condi- 
tions are satisfied, the low-energy physics of the under- 
lying microscopic model can be captured quantitatively 
by the effective field theory. Using the systematic effec- 
tive theory, detailed calculations of a variety of physi- 
cal quantities including 2-loop corrections have been car- 
ried out in Q. Here we only quote the results that arc 
relevant to our study. The aspect ratio of a spatially 
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quadratic space-time box of spatial size L is character- 
ized by I = ((3c/ L) 1 ^ 3 , which distinguishes cubical space- 
time volumes with (3c « L (known as the e-regime in 
QCD) from cylindrical ones with (3c 3> L (the so-called 
S- regime in QCD). In the cubical regime, the volumc- 
and temperature-dependence of the staggered suscepti- 
bility is given by 
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while the uniform susceptibility takes the form 
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In eqs.((5]) and (JB]), the functions (3i(l), (3i(l), and if>(l), 
which only depend on I, are known shape coefficients of 
the space-time box defined in Finally, in the cylin- 
drical regime, when the condition L 2 p a / (f3c 2 ) <C 1 is sat- 
isfied, the volume-dependence of the staggered suscepti- 
bility is given by 
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where a = 0.3103732207 and b = 0.0004304999 [18|. It 
should be noted that Xs given in eq.([7]) is temperature- 
independent. 

Determination of Low-Energy Parameters — In 

order to determine the low-energy parameters A4 S , p s , 
and c for the spin i Heisenberg model on the square lat- 
tice, we have performed large-scale simulations for vari- 
ous inverse temperatures f3 and box sizes L. We deter- 
mine c using the idea proposed in plj : for a fixed box size 
L, we vary (3 until the condition (W t 2 ) = \ {{W?) + (W 2 2 )) 
is satisfied. The spin wave velocity then results as 
c = L//3. Using this method, we obtain c = 1.6586(3)Ja 
(see figure QJ . This value is obtained by performing a 
weighted average over the values of c listed in table I, 
which are extracted in volumes ranging from L = 24a to 
L = 64a. It should be noted that the above value of c is 
consistent with the one quoted in but the statistical 
error is reduced by a factor of 7. In principle, using this 
method one could obtain an even more precise estimate 
of c. After obtaining this very accurate value of c, we 
carry out further large scale simulations in the cubical 
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FIG. 1: The determination of c using the squares of spatial 
and temporal winding numbers at L — 48a. 
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c 


24 


1.6589(6) 


32 


1.6586(5) 


48 


1.6585(5) 


64 


1.6585(5) 



TABLE I: Values of c = L//3 extracted for different lattice 
sizes L/a by tuning f3 such that the average squares of the 
spatial and temporal winding numbers are the same. 



regime with j3c w L. Using c = 1.6586(3)Ja and per- 
forming a combined fit of the Monte Carlo data for \s 
and Xu in the cubical regime to eqs.((S]) and ©, we ar- 
rive at M s = 0.30743(l)/a 2 and p s = 0.18081(11) J with 
X 2 /d.o.f. « 1. Figure [5] illustrates the results of the fit. 
The main contribution to the uncertainties of Ai s and 
p s results from the error of c that enters the fit. Hence, 
with a more precise estimate of c, one could even further 
improve the accuracy of M s and p s . The values we ob- 
tain for p s and c are more accurate than earlier estimates 
of these low-energy parameters. It should be noted that 
the value obtained for M s is consistent with the one of 
[l3T |. and the statistical error is the same in both cases. 

Next we simulate the model in the cylindrical regime 
where the condition (3c 3> L is satisfied. Since a main mo- 
tivation of our study is to clarify the discrepancy between 
the values of M s presented in Q and [l3j|, and the accu- 
racy we must reach is hence below the per mille level, we 
adopt the following strategy. First, we note that at very 
low temperatures Xs becomes temperature- independent. 
In order to avoid underestimating the systematic errors 
in an extrapolation to zero temperature, we simulate at 
sufficiently low temperatures so that Xs becomes inde- 
pendent of (3 within error bars. Second, it should be 
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FIG. 2: Fits of \s and (Wi) (and thus Xu) to their predicted 
behavior in magnon chiral perturbation theory. For better 
visibility, some data used in the fits are omitted in the figure. 



noted that Monte Carlo data for Xs m both the cubical 
Q and the cylindrical regime [3 were used for obtain- 
ing the value of M. s quoted in ['J . Since both p s and 
c obtained in the cubical regime are consistent with the 
corresponding results in the cylindrical regime, one may 
conclude that the overestimation of A4 S = 0.30797(3) / a 2 
presented in Q is due to the cylindrical regime data for 
Xs- In order to minimize statistical correlations between 
Xs and Xu, in our fitting strategy we use only cylindri- 
cal regime data for Xs and cubical regime data for \u- 
Applying these strategies and using c = 1.6586(3) Ja, we 
arrive at M s = 0.30746(4)/a 2 and p s = 0.18081(11)J 
(see figure 3). It should be noted that this value of M Sl 
which we obtain in the cylindrical regime, is consistent 
with the value determined in the cubical regime. It is also 
consistent with the most accurate result for M s that was 
previously obtained [l3| . 

In the cylindrical regime, simulating larger lattices is 
necessary in order to reach the same accuracy for M. s 
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FIG. 3: Fit of Monte Carlo data for Xs m the cylindrical 
regime to their chiral perturbation theory prediction. 



as the one obtained in 13|. This demonstrates the ad- 
vantage of finite-temperature simulations: applying the 
effective field theory predictions to finite-temperature 
data, which can be obtained with a moderate compu- 
tational effort, one achieves a very precise numerical 
value for M. s . Using c = 1. 6586(3) Ja, wc arrive at 
M s = 0.30743(l)/a 2 and p s = 0.18081(10) J from a com- 
bined fit including all available data points. The accuracy 
of these low-energy constants is not improved compared 
to those obtained in the cubical regime alone. This is rea- 
sonable since only a few more data points are included 
in the new fit. 

Finally, we would like to clarify possible reasons for 
the overestimation of M s = 0.30797(3)/a 2 obtained 
in Because of the consistency of both p a and c 

obtained in the cubical and cylindrical regimes d, 0], 
one concludes that the slight overestimation of Ai s 
in [U is due to the cylindrical regime data for \s- 
In particular, in order to employ cq.([7]) to determine 
■Ms, in a Seeley expansion has been performed in 
order to extrapolate the finite-temperature Xs data to 
their corresponding zero-temperature limit. However, 
terminating the Seeley series is a subtle matter. Hence, 
the Seeley extrapolation may lead to an underestimated 
systematic error if the data are outside the window in 
which such an extrapolation is justified. Because of this, 
instead of repeating the analysis performed in Q, we 
adopt another strategy. Specifically, we again simulate 
the model with box sizes L/a = 10, 12, 14, . . . , 20 at 
sufficiently low temperatures so that the Xs data we 
obtain are independent of (3. A combined fit of these 
newly obtained data for Xs at very low temperature 
to eq.Q and the Xu data we obtained earlier in the 
cubical regime to eq.© yields M s = 0.3070(2)/a 2 , 
p s = 0.182(2)J and c = 1.66(1) Ja with x 2 / d -0.f. « 1.3. 
While the obtained value M s = 0.3070(2)/a 2 is slightly 
below A4 S = 0.30743(l)/a 2 , the numerical values of 
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the low-energy parameters that we just obtained are 
indeed consistent with those found in the cubical 
regime calculations [f|. Therefore we conclude that the 
overestimation of M s in Q is most likely due to an 
underestimated systematic error of %s related to the 
termination of the Seeley expansion used in 

Conclusions — In this letter, we have revisited the 
spin h Heisenberg model on the square lattice. In par- 
ticular, we have refined the numerical values of the cor- 
responding low-energy parameters, namely the staggered 
magnetization density M s , the spin stiffness p s , and the 
spin wave velocity c. The spin wave velocity is deter- 
mined to very high accuracy using the squares of spa- 
tial and temporal winding numbers. Remarkably, using 
c = 1.6586(3) Ja together with 2-loop magnon chiral per- 
turbation theory predictions for Xs and \ u in the cu- 
bical regime, we obtained a good fit of more than 150 
data points to two analytic expressions with only two 
unknown parameters. Specifically, from the fit we ob- 
tain M s = 0.30743(l)/a 2 and p s = 0.18081(11)J with 
X 2 /d.o.f. w 1. The precision of M s = 0.30743(l)/a 2 is 
comparable to the one in [l3j which was obtained by an 
unconstrained polynomial fit using up to third or fourth 
powers of l/L. Furthermore, the data used in [l3| were 
obtained at much larger L than those used in our study. 
Thanks to the accurate predictions of magnon chiral per- 
turbation theory, it requires only moderate computing 
resources to reach fraction of a per mille accuracy for 
the low-energy parameters. This is encouraging for QCD 
where such accuracy is mandatory to reach a sufficiently 
precise determination of the fundamental low-energy pa- 
rameters of the strong interaction. If calculations in the 
chiral limit would become feasible, the experience gained 
in the Heisenberg model would suggest that the hyper- 
cubical e-regime would be best suited for extracting the 
low-energy parameters. We have also resolved the puzzle 
of the per mille level discrepancy between the two val- 
ues for M s presented in Q and [u| by re-simulating the 
model in the cylindrical regime. Based on a combined 
fit of Xs in the cylindrical regime and Xu m the cubi- 
cal regime, we arrived at M s = 0.30746(4)/a 2 , which 
is consistent with both the results of [|| and [u|. The 
consistency of the values for M s obtained in the cubical 
and in the cylindrical regime is particularly remarkable 
in view of the increased analytic 2-loop accuracy in the 
cylindrical regime [Hj , which demonstrates the quantita- 
tive correctness of magnon chiral perturbation theory in 



describing the low-energy physics of the underlying mi- 
croscopic model. Finally, we concluded that the small 
discrepancy in the values for M s between Q and [l3| 
may be attributed to the termination of the Seeley ex- 
pansion used in obtaining the former result. 
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